/*
 * get_xplor.c
 *
 *  Created on: 21 avr. 2010
 *      Author: marchi
 */
#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include "get_xplor.h"


/* WARNING:
 * Do _not_ use these routines to calculate the angle between two vectors
 * as acos(cos_angle(u,v)). While it might seem obvious, the acos function
 * is very flat close to -1 and 1, which will lead to accuracy-loss.
 * Instead, use the new gmx_angle() function directly.
 */

int get_xplor(real *** Phi,int nx, int ny, int nz, matrix box, FILE *fp){
	int nx_s,ny_s,nz_s;
	int nx_e,ny_e,nz_e;
	real avg=0.0, avg2=0.0;
	real a=0.0,b=0.0,c=0.0,alpha=0.0, beta=0, gamma=0;

	for(int i=0; i< nx; i++){
		for(int j=0; j< ny; j++){
			for(int k=0; k< nz; k++){
				avg+=Phi[i][j][k];
				avg2+=Phi[i][j][k]*Phi[i][j][k];
			}

		}
	}
	avg/=(real) (nx*ny*nz);
	avg2/=(real) (nx*ny*nz);

	nx_s=0;
	ny_s=0;
	nz_s=0;
	nx_e=nx-1+nx_s;
	ny_e=ny-1+ny_s;
	nz_e=nz-1+nz_s;
	cell_angles(box, &alpha, &beta, &gamma);

	alpha*=RAD2DEG;
	beta*=RAD2DEG;
	gamma*=RAD2DEG;

	a=norm(box[XX]);
	b=norm(box[YY]);
	c=norm(box[ZZ]);
	fprintf(fp,"\n       2 !NTITLE\n") ;
	fprintf(fp," REMARKS Electrostatic Potential from GROMACS\n") ;
	fprintf(fp," REMARKS DATE: 2010-04-21\n") ;
	fprintf(fp,"%8d%8d%8d%8d%8d%8d%8d%8d%8d\n",nx,nx_s,nx_e,ny,ny_s,ny_e,nz,nz_s,nz_e);
	fprintf(fp,"%12.5E%12.5E%12.5E%12.5E%12.5E%12.5E\n",10*a,10*b,10*c,alpha,beta,gamma);
	fprintf(fp, "ZYX\n") ;
	int ja;
	for(int k=0;k<nz;k++){
		fprintf(fp,"%8d\n",k);
		ja=0;
		for(int j=0;j<ny;j++){
			for(int i=0;i<nx;i++){
				ja++;
				fprintf(fp,"%12.5E",Phi[i][j][k]);
				if(!(ja%6)) fprintf(fp,"\n");
			}
		}
		fprintf(fp,"\n");
	}
	fprintf(fp,"%8d\n",-9999);

	real stddev=sqrt(avg2-avg*avg);
	fprintf(fp,"%12.4E %12.4E\n",avg,stddev);

	return 0;
}
